/* This file has the declaration for Tensor4_Riemann*Tensor4 yielding
   a double.  I simplify the expression by removing the identically
   zero components of Riemann. */

/* A(i,j,k,l)*B(i,j,k,l) */

template<class A, class B, char i, char j, char k, char l>
inline double operator*(const Tensor4_Riemann_Expr<A,i,j,k,l> &a,
			const Tensor4_Expr<B,i,j,k,l> &b)
{
  return a(1,0,0,1)*b(1,0,0,1) + a(2,0,0,1)*b(2,0,0,1)
    + a(0,1,0,1)*b(0,1,0,1) + a(2,1,0,1)*b(2,1,0,1)
    + a(0,2,0,1)*b(0,2,0,1) + a(1,2,0,1)*b(1,2,0,1)
    + a(1,0,0,2)*b(1,0,0,2) + a(2,0,0,2)*b(2,0,0,2)
    + a(0,1,0,2)*b(0,1,0,2) + a(2,1,0,2)*b(2,1,0,2)
    + a(0,2,0,2)*b(0,2,0,2) + a(1,2,0,2)*b(1,2,0,2)
    + a(1,0,1,0)*b(1,0,1,0) + a(2,0,1,0)*b(2,0,1,0)
    + a(0,1,1,0)*b(0,1,1,0) + a(2,1,1,0)*b(2,1,1,0)
    + a(0,2,1,0)*b(0,2,1,0) + a(1,2,1,0)*b(1,2,1,0)
    + a(1,0,1,2)*b(1,0,1,2) + a(2,0,1,2)*b(2,0,1,2)
    + a(0,1,1,2)*b(0,1,1,2) + a(2,1,1,2)*b(2,1,1,2)
    + a(0,2,1,2)*b(0,2,1,2) + a(1,2,1,2)*b(1,2,1,2)
    + a(1,0,2,0)*b(1,0,2,0) + a(2,0,2,0)*b(2,0,2,0)
    + a(0,1,2,0)*b(0,1,2,0) + a(2,1,2,0)*b(2,1,2,0)
    + a(0,2,2,0)*b(0,2,2,0) + a(1,2,2,0)*b(1,2,2,0)
    + a(1,0,2,1)*b(1,0,2,1) + a(2,0,2,1)*b(2,0,2,1)
    + a(0,1,2,1)*b(0,1,2,1) + a(2,1,2,1)*b(2,1,2,1)
    + a(0,2,2,1)*b(0,2,2,1) + a(1,2,2,1)*b(1,2,2,1);

//    return a(0,0,0,0)*b(0,0,0,0) + a(1,0,0,0)*b(1,0,0,0) + a(2,0,0,0)*b(2,0,0,0)
//      + a(0,1,0,0)*b(0,1,0,0) + a(1,1,0,0)*b(1,1,0,0) + a(2,1,0,0)*b(2,1,0,0)
//      + a(0,2,0,0)*b(0,2,0,0) + a(1,2,0,0)*b(1,2,0,0) + a(2,2,0,0)*b(2,2,0,0)
//      + a(0,0,0,1)*b(0,0,0,1) + a(1,0,0,1)*b(1,0,0,1) + a(2,0,0,1)*b(2,0,0,1)
//      + a(0,1,0,1)*b(0,1,0,1) + a(1,1,0,1)*b(1,1,0,1) + a(2,1,0,1)*b(2,1,0,1)
//      + a(0,2,0,1)*b(0,2,0,1) + a(1,2,0,1)*b(1,2,0,1) + a(2,2,0,1)*b(2,2,0,1)
//      + a(0,0,0,2)*b(0,0,0,2) + a(1,0,0,2)*b(1,0,0,2) + a(2,0,0,2)*b(2,0,0,2)
//      + a(0,1,0,2)*b(0,1,0,2) + a(1,1,0,2)*b(1,1,0,2) + a(2,1,0,2)*b(2,1,0,2)
//      + a(0,2,0,2)*b(0,2,0,2) + a(1,2,0,2)*b(1,2,0,2) + a(2,2,0,2)*b(2,2,0,2)
//      + a(0,0,1,0)*b(0,0,1,0) + a(1,0,1,0)*b(1,0,1,0) + a(2,0,1,0)*b(2,0,1,0)
//      + a(0,1,1,0)*b(0,1,1,0) + a(1,1,1,0)*b(1,1,1,0) + a(2,1,1,0)*b(2,1,1,0)
//      + a(0,2,1,0)*b(0,2,1,0) + a(1,2,1,0)*b(1,2,1,0) + a(2,2,1,0)*b(2,2,1,0)
//      + a(0,0,1,1)*b(0,0,1,1) + a(1,0,1,1)*b(1,0,1,1) + a(2,0,1,1)*b(2,0,1,1)
//      + a(0,1,1,1)*b(0,1,1,1) + a(1,1,1,1)*b(1,1,1,1) + a(2,1,1,1)*b(2,1,1,1)
//      + a(0,2,1,1)*b(0,2,1,1) + a(1,2,1,1)*b(1,2,1,1) + a(2,2,1,1)*b(2,2,1,1)
//      + a(0,0,1,2)*b(0,0,1,2) + a(1,0,1,2)*b(1,0,1,2) + a(2,0,1,2)*b(2,0,1,2)
//      + a(0,1,1,2)*b(0,1,1,2) + a(1,1,1,2)*b(1,1,1,2) + a(2,1,1,2)*b(2,1,1,2)
//      + a(0,2,1,2)*b(0,2,1,2) + a(1,2,1,2)*b(1,2,1,2) + a(2,2,1,2)*b(2,2,1,2)
//      + a(0,0,2,0)*b(0,0,2,0) + a(1,0,2,0)*b(1,0,2,0) + a(2,0,2,0)*b(2,0,2,0)
//      + a(0,1,2,0)*b(0,1,2,0) + a(1,1,2,0)*b(1,1,2,0) + a(2,1,2,0)*b(2,1,2,0)
//      + a(0,2,2,0)*b(0,2,2,0) + a(1,2,2,0)*b(1,2,2,0) + a(2,2,2,0)*b(2,2,2,0)
//      + a(0,0,2,1)*b(0,0,2,1) + a(1,0,2,1)*b(1,0,2,1) + a(2,0,2,1)*b(2,0,2,1)
//      + a(0,1,2,1)*b(0,1,2,1) + a(1,1,2,1)*b(1,1,2,1) + a(2,1,2,1)*b(2,1,2,1)
//      + a(0,2,2,1)*b(0,2,2,1) + a(1,2,2,1)*b(1,2,2,1) + a(2,2,2,1)*b(2,2,2,1)
//      + a(0,0,2,2)*b(0,0,2,2) + a(1,0,2,2)*b(1,0,2,2) + a(2,0,2,2)*b(2,0,2,2)
//      + a(0,1,2,2)*b(0,1,2,2) + a(1,1,2,2)*b(1,1,2,2) + a(2,1,2,2)*b(2,1,2,2)
//      + a(0,2,2,2)*b(0,2,2,2) + a(1,2,2,2)*b(1,2,2,2) + a(2,2,2,2)*b(2,2,2,2);
}

template<class A, class B, char i, char j, char k, char l>
inline double operator*(const Tensor4_Expr<B,i,j,k,l> &b,
			const Tensor4_Riemann_Expr<A,i,j,k,l> &a)
			
{
  return operator*(a,b);
}
